构建物种丰富度栅格

代码参考：

## 可视化物种丰富度：
# 两种思路：
# 第一种参见：
https://rmacroecology.netlify.app/2018/01/23/a-guide-to-transform-species-shapefiles-into-a-presence-absence-matrix-based-on-a-user-defined-grid-system/
# 第二种参见：
https://rspatial.org/raster/cases/3-speciesdistribution.html

案例

## 学习使用letR绘制物种丰度地图：
library(maps)
map(add=T)

install.packages("letsR")
## 基于R中用户定义的网格将物种分布转换为不存在矩阵的指南
library(letsR)
data("Phyllomedusa")
op <- par(mar = rep(0, 4))  
plot(Phyllomedusa)
## 

## 将物种分布数据转为存在-不存在矩阵；####
# 需要的数据形式是物种名和数据经纬度；
library(spocc)
# Get occurrences for Phyllomedusa
occurrrences <- occ(query = 'Phyllomedusa', from = 'gbif', limit = 5000)
data <- occurrrences$gbif$data$Phyllomedusa
# Remove NAs
remove_na <- is.na(data$longitude) | is.na(data$latitude)
data <- data[!remove_na, ]
xy <- data[, c("longitude", "latitude")] # coordinates
sp_name <- data$name
unique(sp_name)
head(data)
# 现在我们有了坐标和种类名称，我们可以使用该lets.presab.points功能了。
PAM_points <- lets.presab.points(xy, sp_name, xmn = -93, xmx = -29, 
                                 ymn = -57, ymx = 15, res = 1)
plot(PAM_points)
## 基于lets.presab.points函数包含了若干存在-缺失数据代码：
PAM_points$Presence_and_Absence_Matrix
plot(PAM_points$Richness_Raster)
map(add=T)

## 基于构架数据shp来定义物种分布密度 #### 
## 基于物种分布范围的多个shp进行合并后绘制存在-不存在矩阵栅格；
# 注意这里提供的是一个文件夹；
path_Ramphastos <- system.file("extdata", package = "letsR")
rgdaltest<-readOGR('D:/bclang/R-4.0.3/library/letsR/extdata/Ramphastos/Ramphastos_ambiguus_22727999.shp')
PAM_birds <- lets.presab.birds(path_Ramphastos, xmn = -93, xmx = -29, 
                               ymn = -57, ymx = 25)

plot(PAM_birds)

## 设计网格进行投影：####
sp.r <- rasterToPolygons(raster(xmn = -93, xmx = -29,
                                ymn = -57, ymx = 15,
                                resolution = 5))
plot(sp.r)
# Give an ID to the cell
slot(sp.r, "data") <- cbind("ID" = 1:length(sp.r),
                            slot(sp.r, "data"))
plot(sp.r, border = rgb(.5, .5, .5))
map(add = T) ## 这个很好用哎，可以直接进行投影；
## 这里本质上是将182栅格进行叠加：
data(Phyllomedusa)
projection(Phyllomedusa) <- projection(sp.r)
resu <- lets.presab.grid(Phyllomedusa, sp.r, "ID")
rich_plus1 <- rowSums(resu$PAM) + 1
colfunc <- colorRampPalette(c("#fff5f0", "#fb6a4a", "#67000d"))
colors <- c("white", colfunc(max(rich_plus1)))
plot(resu$grid, border = "gray40",
     col = colors[rich_plus1])
map(add = TRUE)

## 将resu$gri由shp转为栅格：
ext <- extent(resu$grid)
xmn <- ext[1]
xmx <- ext[2]
ymn <- ext[3]
ymx <- ext[4]

r <- raster(resolution = 5, xmn = xmn, xmx = xmx,
            ymn = ymn, ymx = ymx)
rich <- rasterize(resu$grid, r, field = "ID")
values(rich) <- rowSums(resu$PAM)
plot(rich)
map(add = TRUE)

r <- raster(ncols=100, nrows=100)
cells <- c(3:5, 210)
r <- rasterFromCells(r, cells)
cbind(1:ncell(r), getValues(r))


r <- raster(ncols=36, nrows=18)
plot(r)
n <- 1000
set.seed(123)
x <- runif(n) * 360 - 180
y <- runif(n) * 180 - 90
dim(xy)
xy <- cbind(x, y)
summary(xy)
# get the (last) indices
r0 <- rasterize(xy, r)
plot(r0)
### 计算物种分布特征 ####
## 计算面积：
# 需要注意的是面积的计算有两种方式，一种是构建最大扩散范围；
# 另外一种是构建点值的buffer；
data("Phyllomedusa")
rangesize <- lets.rangesize(Phyllomedusa,
                            coordinates = "geographic")
rangesize <- rangesize / 1000 # Transform in km2
resu <- lets.maplizer(PAM, rangesize, rownames(rangesize), ras = TRUE)
## 评估纬度和物种分布面积之间的关系：
library(ggplot2)
mpg <- as.data.frame(resu$Matrix)
f <- ggplot(mpg, aes(`Latitude(y)`, Variable_mean))
f + geom_smooth(model = lm) + 
  geom_point(col = rgb(0, 0, 0, .6)) + 
  labs(y = "Range Size") + 
  theme_bw()
## 

## 绘制单个物种的丰富度地图：####
as <- read.csv("D:/xh2/f1_points/xhas.csv") %>% na.omit() %>% data.frame()
as2 = SpatialPoints(cbind(as$longitude, as$latitude), 
                    proj4string = CRS("+init=epsg:4326" ))

## 计算生成矩阵的行和列，使用上面的范围相减；
r <- raster(xmn = 70, xmx = 150,
            ymn = -15, ymx = 45,
            resolution = 5)

sp.r <- rasterToPolygons(raster(xmn = 70, xmx = 150,
                                ymn = -15, ymx = 45,
                                resolution = 5))
##  这样做是对的，只是可能存在可视化错误；
sp <- SpatialPointsDataFrame(as2, as)
rich <- rasterize(sp, r, 'name',function(x, ...) length((x)))
colfunc <- colorRampPalette(c("#fff5f0", "#fb6a4a", "#67000d"))
colors <- c("white", colfunc(max(rich_plus1)))
plot(resu$grid, border = "gray40",
     col = colors[rich_plus1])
plot(rich)
plot(sp.r,add =T,border = "gray40")
points(as2,col="gray70")
map(add =T)

## 用另外一种方法：
r <- raster(xmn = 70, xmx = 150,
            ymn = -15, ymx = 45,
            resolution = 5)
r[] <- 0
tab <- table(cellFromXY(r, as2))     
tab
r[as.numeric(names(tab))] <- tab

plot(sp.r,border = "gray40",add= T)
points(as$longitude,as$latitude,col="grey80")
plot(r,add =T)
map(add = T)